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The accurate representation of geostrophic balance is an essential requirement for numeri- 
cal modelling of geophysical flows. Significant effort is often put into the selection of accu- 
rate or optimal balance representation by the discretisation of the fundamental equations. 
The issue of accurate balance representation is particularly challenging when applying dy- 
namic mesh adaptivity, where there is potential for additional imbalance injection when 
j>. interpolating to new, optimised meshes. 

In the context of shallow-water modelling, we present a new method for preservation of 
Q geostrophic balance when applying dynamic mesh adaptivity. This approach is based upon 

interpolation of the Helmholtz decomposition of the Coriolis acceleration. We apply this in 
combination with a discretisation for which states in geostrophic balance are exactly steady 
solutions of the linearised equations on an /-plane; this method guarantees that a balanced 
and steady flow on a donor mesh remains balanced and steady after interpolation onto an 
arbitrary target mesh, to within machine precision. We further demonstrate the utility of 
this interpolant for states close to geostrophic balance, and show that it prevents pollution 
of the resulting solutions by imbalanced perturbations introduced by the interpolation. 
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1 Introduction 



It has been recognised in previous work on shallow-water modelling that the accurate 
representation of physical balance by the discrete system is of crucial importance. In 



Le Roux et al.| (|1998) a number of shallow- water finite element pairs are compared, and 
it is concluded that the majority of those tested are unable to represent physical balance 
accurately while simultaneously remaining free of spurious numerical modes. In |Cotter 
et al.| ( 2009b[ ) a finite element pair is presented using a piecewise linear discontinuous 
element for velocity and a piecewise quadratic C° continuous element for layer thick- 
ness: the PwgP2 finite element pair. This is shown to be free of spurious pressure modes 
( |Cotter et al.[ [2009b a) and, in |Cotter et al.| ( |2009a| in the context of shallow-water mod- 
elling, is shown to have the property that geostrophically balanced states with a constant 
stream function on the boundary are exactly steady solutions of the discrete linearised 
shallow-water equations on an /-plane. The PidqPi element pair is compared against a 
number of other low order discontinuous methods in |Comblen et al.| (2009), and found 
to be the most accurate choice amongst those tested. This discretisation is extended in 
Cotter et al. ( 2009b|a ) to form a family of related /% ddg^i finite element pairs, all of 
which are free of pressure modes and exhibit this optimal balance property. 



This paper is concerned with the application of dynamic mesh adaptivity to shallow- 
water ocean modelling. Dynamic mesh adaptivity allows the resolution of the computa- 
tional mesh to be varied locally as a simulation develops, in order to resolve dynamically 
important regions of the flow and thereby increase the accuracy per degree of freedom 
of a model. While this has the potential for enabling numerical simulations of other- 
wise inaccessible systems, it presents an additional problem: as well as the possibility 
of imbalance injection in the numerical discretisation of the underlying equations, there 
is also a potential for imbalance injection by the mesh optimisation procedure itself. 



Recently, dynamic mesh adaptive ocean modelling has been proposed by Pain et al. 
(2005); Piggott et al. ( |2008a[ ). This approach utilises unstructured meshes in all three 
dimensions with dynamic mesh adaptivity applied using element-wise topological op- 
erations and nodal perturbations to optimise the mesh according to a metric derived 
from the interpolation error of simulation fields (Pai n et al.[ |2001| |Piggott et al.[ |2006 
Munday et al.[ 12010). A feature of this approach is that, in each mesh adapt, the new 



optimised target mesh has, in general, no relationship to the original donor mesh, other 
than that each is some covering simplex partitioning of the same original domain. This 
generality presents a particularly challenging interpolation problem. 
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The Galerkin projection of fields between arbitrary two and three dimensional meshes 



is described in Farrell et al. (2009); Farrell (2009). This projection is (by definition) 
optimal in the least squares sense, and has the advantage that the integral of the projected 
field is exactly conserved. However, there is no guarantee that the projection injects no 



additional imbalance. In St-Cyr et al. (2008) statically refined mesh and uniform mesh 
simulations are compared in an adaptive mesh refinement (AMR) shallow-water model 
of a geostrophically balanced flow. It is noted that there is an increase in the error in the 
statically refined case, attributed to interpolation errors in the AMR ghost cells which 
"introduce slight disturbance" in this region, with the resulting errors growing faster in 
regions of stronger solution gradients. It is suggested that higher order schemes be used 
to mitigate this. 

In this paper we seek to address the problem of imbalance injection by interpolation 
between arbitrary unstructured meshes. In the context of shallow- water modelling, we 
formulate an interpolant that, for the linearised system on an /-plane, guarantees that 
flows that are initially steady and in geostrophic balance remain steady and in balance 
after interpolation onto an arbitrary target mesh. 



2 Formulation 



In section 2.1 the finite element discretisation of the linearised shallow-water equations 
using the Pidg^2 finite element pair is outlined. A set of properties for geostrophic 



balance preserving interpolants is presented in section 2.2 for which, for the linearised 
system on an /-plane, an initially steady and geostrophically balanced state remains 
steady and balanced after interpolation onto an arbitrary target mesh. An interpolant 



satisfying these properties is given in section 2.3 



2. 1 Discretised shallow -water equations 



The linearised shallow-water equations with free slip boundary conditions are: 



du 

— +fzxu+gVr ] = 0, 



dn 



+ HV • h = 0, 



dt 

u ■ h = on d£l, 



(la) 

(lb) 
(lc) 



where u is the (horizontal) velocity, H is the mean layer thickness, 77 is the deviation 
of the layer thickness from H, f is the Coriolis parameter, g is the gravitational ac- 
celeration and h is a unit normal on the boundary dCl. From this one may define two 
non-dimensional parameters: the Rossby number Ro = U/fD and the Froude number 
Fr = ^JU/gH, where U and D are characteristic flow speeds and spatial scales respec- 
tively. 



Multiplying equations ( |Ta| ) and ([TbJ) by test functions w and respectively, integrating 
over the domain £1, integrating by parts and applying the free slip boundary condition 
yields the weak form: 



J° du 
w — 
a dt 



+ fw ■ (z x u) + gw ■ V77 = Vw, (2a) 

f 4>^--H f V0 • u = V0. (2b) 

Jn dt J a 

Choosing some simplex covering partition of Q. (the mesh), restricting w and u to be 
piecewise linear discontinuous and restricting (p an d rj to be piecewise quadratic C° 
continuous completes the PmoPi spatial discretisation: 



L 

Jn 



Au s 



w-—+fw s -(zxu s ) + gw 6 -Vr] S = Mw s , (3a) 



(V^r -H fv(f> 6 -u 6 = V0 5 , (3b) 

Jn dt Jn 

where £ 6 denotes the finite element approximation for £. Introducing basis function 
expansions of w 6 , u 5 , <f> 6 and j] 6 , this can be re-expressed as: 



—Miii + fLu + gCfj = 0, (4a) 

at 

—M 2 fj - HC T u = 0, (4b) 

dt 



where u and f] are the nodal values for velocity and layer thickness respectively and: 



Mi =diag (M[,M[), 



'o-l^ 



1 



M\, 



M 2 = dmg(M' 2 ,M' 2 ), 



C T = (C\C y ) 



(5) 



M\ and M^ are the velocity space and layer thickness space mass matrices respectively, 
and C T is the discrete divergence matrix: 



(M[) ij = fwj, (M' 2 ).. = ftej, 

(C%= f^j qe{x,y], 



(6) 



where the i/r, and £,• are the Pwg an d Pi elemental basis functions respectively. Choosing 
some time discretisation allows equations (Hal) and (|4b|) to be integrated on a computer. 



2.2 Geostrophic balance preserving interpolants 



Consider interpolation between a donor mesh A and a target mesh B. Let (. . .) A and (. . .) B 
denote "on donor" and "on target" respectively - for example (C T ) A and (C T ) B are the 
divergence matrices, as per Q, assembled on the donor and target meshes respectively. 
Consider an interpolation procedure as follows: 



(1) Perform a Helmholtz decomposition of the Coriolis acceleration F, A = fzxu A on 
the donor mesh: 

M A F A = M A F* A + C A $ A , (7) 

for some scalar potential <D A and discrete divergence free F A . 

(2) Interpolate F A , ® A and f^ from the donor to the target to form F B , O s and fj B , 
using interpolants with the following properties: 



(M A ) l C a O a x h = on dQ => (m b ) 1 C B d> B x h = on dQ, 



F A = => F B = 0, 
o A = gf => o B = ^. 



(8a) 

(8b) 
(8c) 



(3) Recompose F. from O b and F B : 



M B F* B = M B F B - C B d> B . (9) 



The Helmholtz decomposition splits the Coriolis acceleration into a curl-free scalar po- 
tential gradient component and a divergence free residual ( |Weyl[[l9"30" ; Ladyzhenskaya[ 



1969[ ). Only the scalar potential gradient component can be cancelled from equation 



(|4a]) by a layer thickness gradient. For incompressible Navier-Stokes, this scalar po- 
tential gradient component must be exactly cancelled by the pressure gradient, with the 
diagnostic pressure field acting as a Lagrange multiplier via which the incompressibility 
constraint is applied. 



Property ( |8b| ) states the somewhat trivial requirement that zero is preserved by the in- 
terpolant for F. Property ( [8c] ) couples the velocity and layer thickness interpolation, 
and can be achieved if O and fj use the same interpolant and that interpolant is scale 
invariant. Property ([8a]) asserts that the interpolant for O preserves zero tangential gra- 
dients on the domain boundary, and is required in order to avoid generation of grid scale 
boundary Kelvin waves by the interpolation. 

We now proceed to prove that this interpolant, when applied to a /'idg^'2 discretisation 
of the linearised shallow-water equations on an /-plane, guarantees that a steady and 
geostrophically balanced state on the donor mesh results in a state that is steady and 
balanced on the target mesh. By definition, for a geostrophically balanced state on the 
donor: 



fL A u A = -gCV- (10) 



By equation ([7]) F 4 = and <& A = gff, and hence by properties ( [8b] ) and ([8c]) F B = 
and <D B = gf) B . Hence, by equation ([9]), on the target: 



fL B u B = - 8 C B f. (11) 

Also, since F* is perpendicular to u: 



i A ■ h = on dQ 



F/xn 



(Affr'C^D* x n = 0ondQ 
(Mf) _1 C s O B x « = on dQ 
F* x h = on <9Q 



h b • n = on <%1 



by© 

by® 



(12) 



From (Cotter et al. 2009a), if using the PwgPi element pair on an /-plane: 



fLu + gCfj = and u ■ h = on dQ 
dti n dfj 

^¥ = 0and i =a 



(13) 



Hence by (H) , ( [12] ) and ([3\ the solution on the target mesh is geostrophically balanced 
and exactly steady. 



2.3 Implementation and boundary conditions 



The Helmholtz decomposition of the Coriolis acceleration on the donor mesh is equiv- 
alent to the pressure projection method commonly used for incompressible Navier- 
Stokes solvers (Chorin 1967 Temam 1968; Gresho[ 1990). Multiplying equation |7]) 



Tt?A 



by (C y ) (Aff ) and using C 1 F = leads to the elliptic equation: 



(C T ) A (Mfy l C A ® A = -(C T ) A F* 



(14) 



Note that here the consistent mass matrix can be used as the Pidg mass matrix M A is 
block diagonal, and hence the Laplacian matrix (C T ) A (M A )~ l C A is sparse. From this 
O a and F A can be determined. Following interpolation of the Helmholtz decomposition 
F* B can be diagnosed directly from <D B and F B via equation Q. Therefore, the key step 
in forming a geostrophic balance preserving interpolant is to choose interpolants for O, 



f) and F such that the properties ( |8a| ), ( |8b| ) and ( |8c| ) are satisfied 



One simple approach is to apply a Galerkin projection of F from the donor mesh to the 



target mesh, as described in Farrell et al. (2009), and to interpolate <J> and fj using collo- 
cation: evaluation of the donor fields at the nodal coordinates of the target mesh. Since 



Galerkin projection and collocation are linear, properties d8bb and ([8c]) are satisfied. Col 



location also preserves constant boundary values, and hence property ( |8a| ) is satisfied. 
However, collocation (at least for piecewise linear fields) erodes solutions bounds and 



has no optimality properties (Farrell 2009). 



A more accurate approach is to apply a mesh-to-mesh Galerkin projection of O and fj. 
This does not in general satisfy property ( |8a| ), although this issue can be resolved by 
using a further decomposition of O, with an equivalent decomposition of 77 in order to 
satisfy property (f8c|). Assuming Q. is simply connected, O a can be decomposed into: 



A = o£ + ®R> 



(15) 



where O a is equal to some constant c on the boundary, and O^ is some residual. O a can 
be re-expressed: 



A = On + C0i 



(16) 



Here O and Oi satisfy: 



N a O a = 7V a O a , 
JV A Of = 0, 



(17a) 
(17b) 



where N A is the discrete Laplacian matrix N A = (C T ) A (M A ) l C A , N A is N A with a 
Dirichlet boundary condition of zero on dQ. and N A is N A with a Dirichlet boundary 



condition of one on dQ. Minimising O^ 2 subject to ( [15] ) then yields a unique value 
fore: 



(0 A ,O A -O A ) 



L 1 



(of, Of) 



(18) 



L 2 



This choice of c has the property that if O a is constant on dfl, then O^ = 0. Applying 
a Galerkin projection of O^ 
Dirichlet boundary condil 
property ( J8a[ ) is satisfied. 



i A and O^ from the donor mesh to the target mesh, with a 
Dirichlet boundary condition of c on dfl for O a , therefore guarantees that the boundary 



Note that using the mass matrix in place of the Laplacian matrix, /V A = M A , is not 
suitable here, as this results in non-smooth 04 and O^, with strong gradients close 



to the boundary which can generate significant noise in the donor-to-target Galerkin 
projection. 

The full geostrophic balance preserving interpolation procedure is therefore: 

(1) Compute F* A from u A . 

(2) Solve equation ( [14] ) for O 71 and compute F A using equation ([7]). 

(3) Solve equations ( fTTa) and ( fTTb] ), with N A = (C T ) A (M A )~ l C A , for 0>£ and Of, and 
compute 0£ and O^. Perform a similar decomposition for ff" to form rfe and 77^. 

(4) Apply a Galerkin projection from the donor mesh to the target mesh of $£, O^, 77^,, 
77^ and T 7 ' 4 , with Dirichlet boundary conditions for <D^ and 77^. as determined from 



equation ([18]), to form ®£, O^, 77^, 77* and F B . 



(5) Compute O b from <D^ and O^ using equation ( ]T5] ). Similarly compute t/ b from 7)^ 
and 77^. 

(6) Compute F* B using equation ([£[). 

(7) Compute u B from F» B . 



3 Numerical Examples 



In this section several numerical examples of geostrophic balance preservation using 



the interpolation procedure presented above are given. In section 3.1 it is demonstrated 



that the geostrophic balance preserving interpolant ensures that a steady and balanced 



state remains steady and balanced after interpolation. In section 3.2 a state close to 



geostrophic balance is considered, and it is shown that the geostrophic balance preserv- 
ing interpolant avoids imbalance injection. The interpolant is applied to a Kelvin wave 
in section [33] and the accuracy of the interpolant in the L 2 norm is quantified in section 



3.1 Preservation of balance 



The Pidg^2 linearised shallow-water equations d4ab and (4b) on an /-plane were dis 



cretised in time using Crank- Nicolson finite differencing (Crank and Nicolson 1947), 



and the linear systems solved with preconditioned conjugate gradients using the PETSc 
library (Balay et al. , 1997 , 2008, [2009] ) . Further details of the discretisation are given in 
Cotter etal. I ( |2009aD . 



In order to test for imbalance injection by mesh-to-mesh interpolation, two pseudo- 
isotropic circular meshes A and B were generated using Gmsh (Geuzaine and Remacle, 




A B 

Fig. 1. Pseudo-isotropic meshes used to test for imbalance injection by mesh-to-mesh interpola- 
tion. Mesh B has one half the resolution of mesh A. 



2009) and the ani2d mesh optimisation library (VasilevskiiandLipnikov, 1999; Ag ouzal| 
et al. 1999), with mesh A of one half the resolution of mesh B, as shown in figure [T] 
Meshes A and B have 2447 and 557 nodes respectively. Following the balance preser- 
vation test of Le Roux et al. ( 1998); |Cotter et~aL (2009b), the system was initialised 
on mesh A with a Gaussian profile for layer thickness, shown in figure [2} and with a 
velocity field initialised so as to be in discrete geostrophic balance with that profile as 



per equation < |T0] ). The solution was then interpolated backwards and forwards between 
meshes A and B at ten timestep intervals. Since geostrophically balanced states with a 
constant stream function on the boundary are known to be exactly steady when using 



the PivgPi element pair (Cotter et al. , 2009a), any transience observed in the simulation 



is purely due to imbalance injection by the interpolation procedure. 

The model was integrated for 210 timesteps of 6 x 10~ 4 (D/U) at Rossby number 0.06 
and Froude number 0.07 for a total of 20 interpolations between meshes A and B. Three 



interpolants were tested: a first order accurate interpolant as proposed by Grandy ( 1999 1 



for both velocity and layer thickness ("Grandy interpolation"), Galerkin projection for 
velocity and layer thickness, and the geostrophic balance preserving interpolant pre- 
sented in the previous section. Collocation was not tested, as this is not well defined 
at element boundaries and hence is unsuitable for use with the discontinuous velocity 
field. 

The final layer thickness and change in layer thickness from the initial condition are 
shown in figure [3j and the maximum change in layer thickness between each interpo- 
lation is shown in figure |4j Grandy interpolation is observed to inject imbalance every- 
where after each interpolation, resulting in a severe degradation of the simulation fields 



10 




0.0020 0. JO 



Layer Thickness 

0.30 0.40 0.50 0.60 0.70 



Fig. 2. Gaussian profile layer thickness used as an initial condition for the geostrophic balance 
preservation test. 



after just a single interpolation. Galerkin projection is observed to inject imbalance to- 
wards the centre of the domain, near the layer thickness maximum. The resulting gravity 
waves propagate outwards polluting the global solution, and accumulate after every in- 
terpolation. The geostrophic balance preserving interpolant is exactly steady, to within 
machine precision, after every interpolation, with a change in layer thickness between 
interpolations of < 1(T 13 . The residual imbalance between interpolations is attributed to 
double precision round-off error. 

After 20 interpolations the L 2 layer thickness error is 20% (of initial layer thickness 
L 2 norm) for Grandy interpolation, 2.7% for Galerkin projection and 2.0% for the 
geostrophic balance preserving interpolant. While Galerkin projection is optimal in the 
L 2 norm for each interpolation, the imbalance injection and resulting pollution of the 
solution by gravity waves leads to a reduced accuracy in the L 2 norm of the final model 
solution with respect to the geostrophic balance preserving interpolant. 

To further demonstrate geostrophic balance preservation the test was repeated on two 
anisotropic circular meshes C and D generated using Gmsh ( Geuzaine and Remacle[ 
2009) and the ani2d mesh optimisation library (VasilevskiiandLipnikov, 1999; Ag ouzal| 
et al. 1999) with elements stretched in perpendicular directions as shown in figure [5] 
Meshes C and D have 7986 and 7205 nodes respectively, and a maximum element edge 
length ratio of ~ 30. The velocity field was initialised to be in discrete geostrophic 
balance with this layer thickness as before, with interpolations backwards and forwards 
between the two meshes at 10 timestep intervals for 20 interpolations. Simulations were 
conducted using the geostrophic balance preserving interpolant, Galerkin projection, 
and Grandy interpolation, with the change in layer thickness between interpolations 
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Layer Thickness Layer Thickness Error 

0.0037 0.067 0.13 0.19 0.26 0.32 0.39 0.45 0.51 0.58 0.64 -0.079 -0.025 0.029 0.083 0.14 0.J9 0.24 0.30 0.35 0.41 0.46 



II 






Layer Thickness 

0.0019 0.093 0.18 0.27 0.37 0.46 0.55 0.64 0.73 0.82 0.91 - 



III 




Layer Thickness Error 

■0.024 -0.0049 0.014 0.033 0.O52 




Layer Thickness Layer Thickness Error 

0.0 0.094 0.19 0.28 0.37 0.47 0.56 0.66 0.75 0.84 0.94 -0.057 -0.045 -0033 -0.021 -0.0089 0.0030 0.015 0.027 0.039 0.051 0.063 



Fig. 3. The final solution of the geostrophic balance preservation test after 20 repeated interpo- 
lations backwards and forwards between the pseudo-isotropic meshes A and B in figure [T] Left: 
Final layer thickness. Right: Change in layer thickness from the initial condition in figure [2] I: 
Grandy interpolation. II: Galerkin projection. Ill: Helmholtz decomposed geostrophic balance 
preserving interpolation. 

shown in figure [6] When applying the geostrophic balance preserving interpolation the 
maximum change between interpolations was < 10~ 13 as before. 

Finally, the geostrophic balance preservation the test was repeated using the anisotropic 
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Fig. 4. The maximum change in layer thickness between interpolations backwards and forwards 
between the pseudo-isotropic meshes A and B in figure [T] 
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Fig. 5. Anisotropic meshes used to test for imbalance injection by mesh-to-mesh interpolation. 
There is no relationship between meshes C and D, other than that they cover the same domain. 

meshes C and D in figure[5J with a layer thickness initialised to random values in the in- 
terior and a value of zero on the boundary. The velocity field was initialised to be in dis- 
crete geostrophic balance with this layer thickness. The model was integrated as before, 
with interpolations backwards and forwards between two two meshes at 10 timestep 
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Fig. 6. The maximum change in layer thickness between interpolations backwards and forwards 
between the anisotropic meshes C and D in figure [5] 

intervals for 4 interpolations, using the geostrophic balance preserving interpolant. The 
maximum change in layer thickness between interpolations was < 10" 12 , and hence the 
solution was observed to be steady to within double precision round-off error. 



3.2 Nearly balanced states 



The Gaussian layer thickness profile in figure [2]had a perturbation applied of the form: 



Ai/ = ^Xfj, 



(19) 



where X is some point-wise random value in the range {0 - 1}. This perturbation was 
smoothed using a Helmholtz smoother with a characteristic length scale of D/8 to pro- 
duce the layer thickness perturbation shown in figure|7J The velocity field was initialised 
to be in discrete geostrophic balance with the unperturbed layer thickness, thereby gen- 
erating a nearly balanced state. 

The system was integrated as before, with interpolations backwards and forwards be- 
tween the pseudo-isotropic meshes A and B in figure[T]at 10 timestep intervals. One can 
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Layer Thickness Perturbation 

0.0000 0.0004 0.0008 0.0013 0.0017 0.0021 0.0025 0.0029 



Fig. 7. Perturbation applied to the layer thickness in figure[2]to test for imbalance injection by 
interpolation of a nearly balanced state. 

define an "imbalanced layer thickness": 



fjimbal -=n- ~®, (20) 

8 

where <D is the scalar potential computed from the Helmholtz decomposition of the Cori- 
olis acceleration. The final imbalanced layer thickness is shown for Galerkin projection 
and the geostrophic balance preserving interpolant in figure [8] When using Galerkin 
projection imbalance is observed to be injected near the layer thickness maximum. This 
additional imbalance dominates over the original layer thickness perturbation after 20 
interpolations. When using the geostrophic balance preserving interpolant propagation 
of the original layer thickness perturbation is observed, with no significant pollution 
introduced by the interpolation. 

Defining a "balanced velocity" Ubai where: 

fl&bal '■= -gCfj, (21) 

and an "imbalanced velocity" u im b a i- 

Uimbal '■= U — Ubah (22) 

allows one to compute an imbalanced kinetic energy: 

15 





Imbalanced Layer Thickness Imbalanced Layer Thickness 

-0.625 -0.677 -0.529 -0.381 -0.233 -0.0844 0.0636 0.212 0.360 0.508 0.656 -0.158 -0.137 -0.116 -0.0943 -0.0730 -0.OS17 -0.0304 -0.00909 0.0122 0.0335 0.0548 



A 



B 



Fig. 8. The final imbalanced layer thickness, computed as the difference between the simulation 
layer thickness and the scaled scalar potential from the Helmholtz decomposition of the Coriolis 
acceleration. A: Galerkin projection. B: Helmholtz decomposed geostrophic balance preserving 
interpolation. 



-* imbal — 2 \\Mimbal\\[2 



a + -L^cn 
f 



(23) 



The imbalanced kinetic energies when using Galerkin projection and the 
geostrophic balance preserving interpolant are shown in figure [9] When using Galerkin 
projection the imbalanced kinetic energy is observed to increase by up to a factor of 70 
in an interpolation, with the imbalanced kinetic energy peaking at 150 times its initial 
value. When using the geostrophic balance preserving interpolant the imbalanced ki- 
netic energy is observed to increase by at most a factor 1 .02 in an interpolation, and the 
imbalanced kinetic energy never exceeds its initial value. 
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Fig. 9. The imbalanced kinetic energy, normalised by the initial imbalanced kinetic energy, for 
the nearly balanced interpolation test. When using Galerkin projection (upper line) large in- 
creases in the imbalanced kinetic energy are observed after interpolation, with these increases 
dominating over the original imbalanced kinetic energy. When using the geostrophic balance 
preserving interpolant (lower line) imbalanced kinetic energy injection is significantly reduced. 

3.3 Kelvin wave 



The interpolant was tested for a Kelvin wave, configured as in in Ham et al. (20051; 



Cotter et al. (2009b I with an initial condition: 



r] (r, 6) = exp I ) cos 6, 

\ Ro 

1 ( r ~ r o\ 

u e (r, 0) = — exp I — — ] cos G, 

u r = 0, 



(24a) 

(24b) 
(24c) 



for Ro = 10 and Fr = 1 in a circular domain of radius r . The Kelvin wave is geostroph- 
ically balanced in the direction normal to the boundary and imbalanced in the tangential 
direction. The model was integrated with a timestep of In x 10~ 4 (D/U) for a total sim- 
ulated time of 2n(D/U), corresponding to the time taken for a single Kelvin wave to 
perform a circuit of the domain in the limit of large r . Two meshes of quasi-uniform 
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Fig. 10. Initial layer thickness used for the Kelvin wave test at Ro = 10, Fr - 1. 



resolution with 1473 and 1461 nodes respectively were created using gmsh ( jGeuzaine] 
and Remacle , 2009 ) and the ani2d mesh optimisation library ( Vas ilevskii and Lipnikov[ 



1999 Agouzal et al. 1999), and the solution interpolated backwards and forwards be- 



tween these meshes at 10 timestep intervals. 



The initial layer thickness is shown in figure 10 and the Helmholtz decomposition of 



the initial Coriolis acceleration in figure 1 1 The final solutions when using Galerkin 



projection and the geostrophic balance preserving interpolant are shown in figure 12 



Relatively little difference is observed in the final layer thickness field between these 
simulations. However, when using Galerkin projection, noise is observed in the velocity 
divergence field, originating at the boundary. This noise is significantly reduced when 
using the geostrophic balance preserving interpolant. 

A discretisation of the linearised shallow-water equations conserves energy if the layer 
thickness gradient matrix is, after multiplication by some diagonal matrix, equal to the 



transpose of the velocity divergence matrix (Ham et al.. 2007), and if the implicit mid 



point rule is used for timestepping (Leimkuhler and Reich 2004). Hence the P\t>gP2 
spatial discretisation of the linearised shallow-water equations as presented here con- 
serves the total energy. The kinetic, potential, and total energy of the system when using 
Galerkin projection, the geostrophic balance preserving interpolant, and when using a 



single fixed computational mesh, are shown in figure 13 The fixed mesh simulation is 
observed to conserve the total energy to within one part in 10 4 , with the relatively high 
error attributed to the tolerances used for the linear solvers. The use of direct solvers, 
combined with more precision robust calculation of the energy diagnostics, is expected 
to decrease this error. When interpolating between meshes using Galerkin projection 
a systematic dissipation of both kinetic and potential energy is observed, leading to a 
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Fig. 11. Helmholtz decomposition of initial Coriolis acceleration for the Kelvin wave test at 
Ro = 10, Fr - 1. A: Coriolis acceleration, F*. B: Non-divergent residual, F. C: Scalar potential 
with a constant boundary value, O c . D: Scalar potential residual, ®#. 

decrease in the total system energy of 1.3% after 1000 interpolations, at the end of the 
simulation. When using the geostrophic balance preserving interpolant a slight increase 
the potential energy is observed, leading to an increase in the total system energy of 
0.11% after 1000 interpolations. While the geostrophic balance preserving interpolant 
is not energy conserving, the change in system energy is, for this test, more than an 
order of magnitude smaller than that observed when applying Galerkin projection. 

In further testing it was found that highly anisotropic elements intersecting the domain 
boundary led to very poor results when using the geostrophic balance preserving in- 
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Fig. 12. The final solution of the Kelvin wave test at Ro = 10, Fr - 1 Left: Final layer thickness. 
Right: Final velocity divergence, M~ l C T u. A: Galerkin projection. B: Helmholtz decomposed 
geostrophic balance preserving interpolation. 

terpolant. This is likely due to significant interpolation errors in the projection of O in 
this region, possibly as a result of the Dirichlet boundary condition for <D C in the mesh- 
to-mesh Galerkin projection, which pollutes the interpolated Coriolis acceleration. This 
problem was solved by imposing constraints on the maximum element size for elements 
directly on the boundary. 



3.4 Accuracy 



Since the geostrophic balance preserving interpolant is composed of a Galerkin dis- 
cretisation of the Helmholtz decomposition of the Coriolis acceleration followed by a 
donor-to-target Galerkin projection of the decomposition, when using the P\dgP2 ele- 
ment pair the interpolant is expected to be second order accurate for velocity and third 
order accurate for layer thickness. 
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Fig. 13. The energy of the Kelvin wave test when interpolating between meshes using Galerkin 
projection and the geostrophic balance preserving interpolant, and when using a single fixed 
computational mesh. A: Kinetic energy. B: Potential energy. C: Total energy. 

A series of structured triangular mesh pairs for a 2D unit square -0.5 < x < 0.5 and 
-0.5 < y < 0.5 were generated with resolutions in the x- and y-directions as given in 
table [TJ A layer thickness of the form: 



7] = sin(2.57u;)sin(2.57r;y), 



(25) 
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Donor mesh resolution (x x y) 


Target mesh resolution(x x y) 


24x26 


26x24 


32x35 


35x32 


48x52 


52x48 


64x70 


70x64 


96 x 105 


105 x 96 


128 x 130 


130x128 


192x215 


215 x 196 



Table 1 

Mesh resolutions used for the geostrophic balance preserving interpolant convergence test. A 
mesh resolution of N x M denotes a division in the x-direction into N sections, a division in the 
y-direction into M sections, and the division of each resulting quadrilateral into two triangles. 



and a velocity of the form: 



u = -z x V?7 + sin (0.5x)jt + sin(0.5x)j, 



(26) 



were interpolated between the meshes in each pair using the geostrophic balance pre- 
serving interpolant. The first term in p.4| ) corresponds to a flow that is, for / = -1, g = 1 



and H = 1, in geostrophic balance with the layer thickness ( |3.4[ ). The remaining terms 
correspond to a flow that cannot be balanced by any layer thickness. In order to test 
for additional error introduced by the scalar potential O and layer thickness decomposi- 



tion, as per equation ( [15] ), tests were conducted for a doubly periodic and for a bounded 
domain. The L 2 errors \\fj B - fj A \\ L 2, \\u x ,b ~ u x a\\ L 2 an d \\i*y,B ~ Uyj\\ T 2> were computed 
explicitly via supermesh construction, as described in Farrell (2009). For comparison 
the fields were also projected using Galerkin projection, giving a measure of the qual- 
ity of the geostrophic balance preserving interpolant relative to the projection that is 
optimal in the L 2 norm. 



The resulting errors are shown for the doubly periodic domain in figure 14 and for 



the bounded domain in figure 15 The geostrophic balance preserving interpolant is 



observed to be second order accurate for velocity and third order accurate for layer 
thickness, as expected. For the doubly periodic domain the average L 2 norm error for 
the geostrophic balance preserving interpolant is observed to be 1 .37 times the optimal 
value for velocity, and (since no layer thickness decomposition is applied in this case) 
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Fig. 14. Convergence test for the geostrophic balance preserving interpolant in a doubly periodic 
domain, with Galerkin projection for comparison. Left: L 2 error in the layer thickness. Right: L 2 
error in the x-component of velocity. The error in the y-component of velocity is similar. 
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Fig. 15. Convergence test for the geostrophic balance preserving interpolant in a bounded do- 
main, with Galerkin projection for comparison. Left: L 2 error in the layer thickness. Right: L2 
error in the x-component of velocity. The error in the y-component of velocity is similar. 



optimal for layer thickness. For the bounded domain the error in velocity is not signifi- 
cantly changed, and the error in layer thickness is increased to 1 .005 times the optimal 
value, indicating that the decomposition of the scalar potential <D A and layer thickness f) A 
introduces only a small additional error. For comparison, in Farrell ( |2009 ) collocation 
is found to give, for a field sin x + cos x, an L 2 error that is ~ 2 - 2.5 times the optimal 
value for piecewise linear elements, and ~ 1,1 times the optimal value for piecewise 
quadratic elements. 
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4 Conclusions 



We have presented an interpolation method that, when applied to the PydgPi discretisa- 
tion of the linearised shallow-water equations on an /-plane, guarantees that steady and 
geostrophically balanced states on the donor mesh remain steady and geostrophically 
balanced after interpolation onto an arbitrary target mesh. We have stress tested this bal- 
ance preserving property with highly anisotropic meshes and randomly initialised bal- 
anced states (constrained to satisfy appropriate boundary conditions). We have further 
demonstrated the utility of this interpolant for nearly balanced dynamics, and quantified 
its accuracy in the L 2 norm. 



A shortcoming of this approach, at least in the form presented, is that is does not con- 
serve energy. The Helmholtz decomposed interpolation of Coriolis acceleration does 
not conserve kinetic energy or potential energy. Despite this, the change in energy when 
using the geostrophic balance preserving interpolant was found to be more than an or- 
der of magnitude smaller than the energy dissipation when using Galerkin projection. 
In addition to this, the interpolant does not locally conserve potential vorticity. Geo- 
physical flows are only in geostrophic balance to leading order, and a lack of potential 
vorticity conservation could, when non-linear advection is included in the fundamental 
equations, lead to higher order balance loss. Potential vorticity decompositions could be 
considered where such a conservation is desired (Sta quet and Riley[|1989^|Holopainen| 
and Kaurola[ 1991, Mclntyre and Norton 2000| ), although the benefit of such an ap- 
proach, bearing in mind that discretisations of the non-linear shallow-water equations 
are not generally potential vorticity conserving, may be somewhat limited compared to 
the benefit of leading order balance preservation. 



The method has a natural extension to Navier-Stokes. For incompressible Navier-Stokes 
any forcing that can be represented as the gradient of a scalar field must be filtered by 
the pressure gradient, and hence the interpolation of the Helmholtz decomposition of the 
Coriolis acceleration is a balance preserving interpolant. Future work will concentrate 
on the implementation of geostrophic balance preserving interpolation as part of the Im- 
perial College Ocean Model - an unstructured dynamic mesh adaptive ocean model. In 
particular, we will investigate accurate preservation of geostrophic balance when using 
velocity-pressure element pairs that do not satisfy optimal balance properties, and the 
integration of methods used for accurate balance representation for such element pairs 
( |Ford et al.[ |2004a|b[ |Piggott et si] [2006] |2008b[ pang et al.[ [2009]) into a balance pre- 
serving interpolant. We will test how these can be used to propagate accurate balance 
representation through arbitrary mesh adapts for meshes that are fully unstructured in 
all three dimensions. 
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